home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / rng / random.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-12-09  |  15.6 KB  |  657 lines
  1. /* rng/random.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_rng.h>
  23.  
  24. /* This file provides support for random() generators. There are three
  25.    versions in widespread use today,
  26.  
  27.    - The original BSD version, e.g. on SunOS 4.1 and FreeBSD.
  28.  
  29.    - The Linux libc5 version, which is differs from the BSD version in
  30.    its seeding procedure, possibly due to the introduction of a typo
  31.    in the multiplier.
  32.  
  33.    - The GNU glibc2 version, which has a new (and better) seeding
  34.    procedure.
  35.  
  36.    They all produce different numbers, due to the different seeding
  37.    algorithms, but the algorithm for the generator is the same in each
  38.    case.
  39.   
  40.  */
  41.  
  42. static inline long int random_get (int * i, int * j, int n, long int * x);
  43.  
  44. static inline unsigned long int random8_get (void *vstate);
  45. static inline unsigned long int random32_get (void *vstate);
  46. static inline unsigned long int random64_get (void *vstate);
  47. static inline unsigned long int random128_get (void *vstate);
  48. static inline unsigned long int random256_get (void *vstate);
  49.  
  50. static double random8_get_double (void *vstate);
  51. static double random32_get_double (void *vstate);
  52. static double random64_get_double (void *vstate);
  53. static double random128_get_double (void *vstate);
  54. static double random256_get_double (void *vstate);
  55.  
  56. static void random8_glibc2_set (void *state, unsigned long int s);
  57. static void random32_glibc2_set (void *state, unsigned long int s);
  58. static void random64_glibc2_set (void *state, unsigned long int s);
  59. static void random128_glibc2_set (void *state, unsigned long int s);
  60. static void random256_glibc2_set (void *state, unsigned long int s);
  61.  
  62. static void random8_libc5_set (void *state, unsigned long int s);
  63. static void random32_libc5_set (void *state, unsigned long int s);
  64. static void random64_libc5_set (void *state, unsigned long int s);
  65. static void random128_libc5_set (void *state, unsigned long int s);
  66. static void random256_libc5_set (void *state, unsigned long int s);
  67.  
  68. static void random8_bsd_set (void *state, unsigned long int s);
  69. static void random32_bsd_set (void *state, unsigned long int s);
  70. static void random64_bsd_set (void *state, unsigned long int s);
  71. static void random128_bsd_set (void *state, unsigned long int s);
  72. static void random256_bsd_set (void *state, unsigned long int s);
  73.  
  74. static void bsd_initialize (long int * x, int n, unsigned long int s);
  75. static void libc5_initialize (long int * x, int n, unsigned long int s);
  76. static void glibc2_initialize (long int * x, int n, unsigned long int s);
  77.  
  78. typedef struct
  79.   {
  80.     long int x;
  81.   }
  82. random8_state_t;
  83.  
  84. typedef struct
  85.   {
  86.     int i, j;
  87.     long int x[7];
  88.   }
  89. random32_state_t;
  90.  
  91. typedef struct
  92.   {
  93.     int i, j;
  94.     long int x[15];
  95.   }
  96. random64_state_t;
  97.  
  98. typedef struct
  99.   {
  100.     int i, j;
  101.     long int x[31];
  102.   }
  103. random128_state_t;
  104.  
  105. typedef struct
  106.   {
  107.     int i, j;
  108.     long int x[63];
  109.   }
  110. random256_state_t;
  111.  
  112. static inline unsigned long int
  113. random8_get (void *vstate)
  114. {
  115.   random8_state_t *state = (random8_state_t *) vstate;
  116.  
  117.   state->x = (1103515245 * state->x + 12345) & 0x7fffffffUL;
  118.   return state->x;
  119. }
  120.  
  121. static inline long int
  122. random_get (int * i, int * j, int n, long int * x)
  123. {
  124.   long int k ;
  125.  
  126.   x[*i] += x[*j] ;
  127.   k = (x[*i] >> 1) & 0x7FFFFFFF ;
  128.   
  129.   (*i)++ ;
  130.   if (*i == n)
  131.     *i = 0 ;
  132.   
  133.   (*j)++ ;
  134.   if (*j == n)
  135.     *j = 0 ;
  136.  
  137.   return k ;
  138. }
  139.  
  140. static inline unsigned long int
  141. random32_get (void *vstate)
  142. {
  143.   random32_state_t *state = (random32_state_t *) vstate;
  144.   unsigned long int k = random_get (&state->i, &state->j, 7, state->x) ; 
  145.   return k ;
  146. }
  147.  
  148. static inline unsigned long int
  149. random64_get (void *vstate)
  150. {
  151.   random64_state_t *state = (random64_state_t *) vstate;
  152.   long int k = random_get (&state->i, &state->j, 15, state->x) ; 
  153.   return k ;
  154. }
  155.  
  156. static inline unsigned long int
  157. random128_get (void *vstate)
  158. {
  159.   random128_state_t *state = (random128_state_t *) vstate;
  160.   unsigned long int k = random_get (&state->i, &state->j, 31, state->x) ; 
  161.   return k ;
  162. }
  163.  
  164. static inline unsigned long int
  165. random256_get (void *vstate)
  166. {
  167.   random256_state_t *state = (random256_state_t *) vstate;
  168.   long int k = random_get (&state->i, &state->j, 63, state->x) ; 
  169.   return k ;
  170. }
  171.  
  172. static double
  173. random8_get_double (void *vstate)
  174. {
  175.   return random8_get (vstate) / 2147483648.0 ;
  176. }
  177.  
  178. static double
  179. random32_get_double (void *vstate)
  180. {
  181.   return random32_get (vstate) / 2147483648.0 ;
  182. }
  183.  
  184. static double
  185. random64_get_double (void *vstate)
  186. {
  187.   return random64_get (vstate) / 2147483648.0 ;
  188. }
  189.  
  190. static double
  191. random128_get_double (void *vstate)
  192. {
  193.   return random128_get (vstate) / 2147483648.0 ;
  194. }
  195.  
  196. static double
  197. random256_get_double (void *vstate)
  198. {
  199.   return random256_get (vstate) / 2147483648.0 ;
  200. }
  201.  
  202. static void
  203. random8_bsd_set (void *vstate, unsigned long int s)
  204. {
  205.   random8_state_t *state = (random8_state_t *) vstate;
  206.   
  207.   if (s == 0) 
  208.     s = 1;
  209.  
  210.   state->x = s;
  211. }
  212.  
  213. static void
  214. random32_bsd_set (void *vstate, unsigned long int s)
  215. {
  216.   random32_state_t *state = (random32_state_t *) vstate;
  217.   int i;
  218.  
  219.   bsd_initialize (state->x, 7, s) ;
  220.  
  221.   state->i = 3;
  222.   state->j = 0;
  223.   
  224.   for (i = 0 ; i < 10 * 7 ; i++)
  225.     random32_get (state) ; 
  226. }
  227.  
  228. static void
  229. random64_bsd_set (void *vstate, unsigned long int s)
  230. {
  231.   random64_state_t *state = (random64_state_t *) vstate;
  232.   int i;
  233.  
  234.   bsd_initialize (state->x, 15, s) ;
  235.  
  236.   state->i = 1;
  237.   state->j = 0;
  238.   
  239.   for (i = 0 ; i < 10 * 15 ; i++)
  240.     random64_get (state) ; 
  241. }
  242.  
  243. static void
  244. random128_bsd_set (void *vstate, unsigned long int s)
  245. {
  246.   random128_state_t *state = (random128_state_t *) vstate;
  247.   int i;
  248.  
  249.   bsd_initialize (state->x, 31, s) ;
  250.  
  251.   state->i = 3;
  252.   state->j = 0;
  253.   
  254.   for (i = 0 ; i < 10 * 31 ; i++)
  255.     random128_get (state) ; 
  256. }
  257.  
  258. static void
  259. random256_bsd_set (void *vstate, unsigned long int s)
  260. {
  261.   random256_state_t *state = (random256_state_t *) vstate;
  262.   int i;
  263.  
  264.   bsd_initialize (state->x, 63, s) ;
  265.  
  266.   state->i = 1;
  267.   state->j = 0;
  268.   
  269.   for (i = 0 ; i < 10 * 63 ; i++)
  270.     random256_get (state) ; 
  271. }
  272.  
  273. static void 
  274. bsd_initialize (long int * x, int n, unsigned long int s)
  275. {
  276.   int i; 
  277.  
  278.   if (s == 0)
  279.     s = 1 ;
  280.  
  281.   x[0] = s;
  282.  
  283.   for (i = 1 ; i < n ; i++)
  284.     x[i] = 1103515245 * x[i-1] + 12345 ;
  285. }
  286.  
  287. static void 
  288. libc5_initialize (long int * x, int n, unsigned long int s)
  289. {
  290.   int i; 
  291.  
  292.   if (s == 0)
  293.     s = 1 ;
  294.  
  295.   x[0] = s;
  296.  
  297.   for (i = 1 ; i < n ; i++)
  298.     x[i] = 1103515145 * x[i-1] + 12345 ;
  299. }
  300.  
  301. static void 
  302. glibc2_initialize (long int * x, int n, unsigned long int s)
  303. {
  304.   int i; 
  305.  
  306.   if (s == 0)
  307.     s = 1 ;
  308.  
  309.   x[0] = s;
  310.  
  311.   for (i = 1 ; i < n ; i++)
  312.     {
  313.       const long int h = s / 127773;
  314.       const long int t = 16807 * (s - h * 127773) - h * 2836;
  315.       if (t < 0)
  316.     {
  317.       s = t + 2147483647 ;
  318.     }
  319.       else
  320.     {
  321.       s = t ;
  322.     }
  323.  
  324.     x[i] = s ;
  325.     }
  326. }
  327.  
  328. static void
  329. random8_glibc2_set (void *vstate, unsigned long int s)
  330. {
  331.   random8_state_t *state = (random8_state_t *) vstate;
  332.   
  333.   if (s == 0) 
  334.     s = 1;
  335.  
  336.   state->x = s;
  337. }
  338.  
  339. static void
  340. random32_glibc2_set (void *vstate, unsigned long int s)
  341. {
  342.   random32_state_t *state = (random32_state_t *) vstate;
  343.   int i;
  344.  
  345.   glibc2_initialize (state->x, 7, s) ;
  346.  
  347.   state->i = 3;
  348.   state->j = 0;
  349.   
  350.   for (i = 0 ; i < 10 * 7 ; i++)
  351.     random32_get (state) ; 
  352. }
  353.  
  354. static void
  355. random64_glibc2_set (void *vstate, unsigned long int s)
  356. {
  357.   random64_state_t *state = (random64_state_t *) vstate;
  358.   int i;
  359.  
  360.   glibc2_initialize (state->x, 15, s) ;
  361.  
  362.   state->i = 1;
  363.   state->j = 0;
  364.   
  365.   for (i = 0 ; i < 10 * 15 ; i++)
  366.     random64_get (state) ; 
  367. }
  368.  
  369. static void
  370. random128_glibc2_set (void *vstate, unsigned long int s)
  371. {
  372.   random128_state_t *state = (random128_state_t *) vstate;
  373.   int i;
  374.  
  375.   glibc2_initialize (state->x, 31, s) ;
  376.  
  377.   state->i = 3;
  378.   state->j = 0;
  379.   
  380.   for (i = 0 ; i < 10 * 31 ; i++)
  381.     random128_get (state) ; 
  382. }
  383.  
  384. static void
  385. random256_glibc2_set (void *vstate, unsigned long int s)
  386. {
  387.   random256_state_t *state = (random256_state_t *) vstate;
  388.   int i;
  389.  
  390.   glibc2_initialize (state->x, 63, s) ;
  391.  
  392.   state->i = 1;
  393.   state->j = 0;
  394.   
  395.   for (i = 0 ; i < 10 * 63 ; i++)
  396.     random256_get (state) ; 
  397. }
  398.  
  399.  
  400. static void
  401. random8_libc5_set (void *vstate, unsigned long int s)
  402. {
  403.   random8_state_t *state = (random8_state_t *) vstate;
  404.   
  405.   if (s == 0) 
  406.     s = 1;
  407.  
  408.   state->x = s;
  409. }
  410.  
  411. static void
  412. random32_libc5_set (void *vstate, unsigned long int s)
  413. {
  414.   random32_state_t *state = (random32_state_t *) vstate;
  415.   int i;
  416.  
  417.   libc5_initialize (state->x, 7, s) ;
  418.  
  419.   state->i = 3;
  420.   state->j = 0;
  421.   
  422.   for (i = 0 ; i < 10 * 7 ; i++)
  423.     random32_get (state) ; 
  424. }
  425.  
  426. static void
  427. random64_libc5_set (void *vstate, unsigned long int s)
  428. {
  429.   random64_state_t *state = (random64_state_t *) vstate;
  430.   int i;
  431.  
  432.   libc5_initialize (state->x, 15, s) ;
  433.  
  434.   state->i = 1;
  435.   state->j = 0;
  436.   
  437.   for (i = 0 ; i < 10 * 15 ; i++)
  438.     random64_get (state) ; 
  439. }
  440.  
  441. static void
  442. random128_libc5_set (void *vstate, unsigned long int s)
  443. {
  444.   random128_state_t *state = (random128_state_t *) vstate;
  445.   int i;
  446.  
  447.   libc5_initialize (state->x, 31, s) ;
  448.  
  449.   state->i = 3;
  450.   state->j = 0;
  451.   
  452.   for (i = 0 ; i < 10 * 31 ; i++)
  453.     random128_get (state) ; 
  454. }
  455.  
  456. static void
  457. random256_libc5_set (void *vstate, unsigned long int s)
  458. {
  459.   random256_state_t *state = (random256_state_t *) vstate;
  460.   int i;
  461.  
  462.   libc5_initialize (state->x, 63, s) ;
  463.  
  464.   state->i = 1;
  465.   state->j = 0;
  466.   
  467.   for (i = 0 ; i < 10 * 63 ; i++)
  468.     random256_get (state) ; 
  469. }
  470.  
  471. static const gsl_rng_type random_glibc2_type =
  472. {"random-glibc2",            /* name */
  473.  0x7fffffffUL,            /* RAND_MAX */
  474.  0,                /* RAND_MIN */
  475.  sizeof (random128_state_t),
  476.  &random128_glibc2_set,
  477.  &random128_get,
  478.  &random128_get_double};
  479.  
  480. static const gsl_rng_type random8_glibc2_type =
  481. {"random8-glibc2",            /* name */
  482.  0x7fffffffUL,            /* RAND_MAX */
  483.  0,                /* RAND_MIN */
  484.  sizeof (random8_state_t),
  485.  &random8_glibc2_set,
  486.  &random8_get,
  487.  &random8_get_double};
  488.  
  489. static const gsl_rng_type random32_glibc2_type =
  490. {"random32-glibc2",            /* name */
  491.  0x7fffffffUL,            /* RAND_MAX */
  492.  0,                /* RAND_MIN */
  493.  sizeof (random32_state_t),
  494.  &random32_glibc2_set,
  495.  &random32_get,
  496.  &random32_get_double};
  497.  
  498. static const gsl_rng_type random64_glibc2_type =
  499. {"random64-glibc2",            /* name */
  500.  0x7fffffffUL,            /* RAND_MAX */
  501.  0,                /* RAND_MIN */
  502.  sizeof (random64_state_t),
  503.  &random64_glibc2_set,
  504.  &random64_get,
  505.  &random64_get_double};
  506.  
  507. static const gsl_rng_type random128_glibc2_type =
  508. {"random128-glibc2",            /* name */
  509.  0x7fffffffUL,            /* RAND_MAX */
  510.  0,                /* RAND_MIN */
  511.  sizeof (random128_state_t),
  512.  &random128_glibc2_set,
  513.  &random128_get,
  514.  &random128_get_double};
  515.  
  516. static const gsl_rng_type random256_glibc2_type =
  517. {"random256-glibc2",            /* name */
  518.  0x7fffffffUL,            /* RAND_MAX */
  519.  0,                /* RAND_MIN */
  520.  sizeof (random256_state_t),
  521.  &random256_glibc2_set,
  522.  &random256_get,
  523.  &random256_get_double};
  524.  
  525. static const gsl_rng_type random_libc5_type =
  526. {"random-libc5",            /* name */
  527.  0x7fffffffUL,            /* RAND_MAX */
  528.  0,                /* RAND_MIN */
  529.  sizeof (random128_state_t),
  530.  &random128_libc5_set,
  531.  &random128_get,
  532.  &random128_get_double};
  533.  
  534. static const gsl_rng_type random8_libc5_type =
  535. {"random8-libc5",            /* name */
  536.  0x7fffffffUL,            /* RAND_MAX */
  537.  0,                /* RAND_MIN */
  538.  sizeof (random8_state_t),
  539.  &random8_libc5_set,
  540.  &random8_get,
  541.  &random8_get_double};
  542.  
  543. static const gsl_rng_type random32_libc5_type =
  544. {"random32-libc5",            /* name */
  545.  0x7fffffffUL,            /* RAND_MAX */
  546.  0,                /* RAND_MIN */
  547.  sizeof (random32_state_t),
  548.  &random32_libc5_set,
  549.  &random32_get,
  550.  &random32_get_double};
  551.  
  552. static const gsl_rng_type random64_libc5_type =
  553. {"random64-libc5",            /* name */
  554.  0x7fffffffUL,            /* RAND_MAX */
  555.  0,                /* RAND_MIN */
  556.  sizeof (random64_state_t),
  557.  &random64_libc5_set,
  558.  &random64_get,
  559.  &random64_get_double};
  560.  
  561. static const gsl_rng_type random128_libc5_type =
  562. {"random128-libc5",            /* name */
  563.  0x7fffffffUL,            /* RAND_MAX */
  564.  0,                /* RAND_MIN */
  565.  sizeof (random128_state_t),
  566.  &random128_libc5_set,
  567.  &random128_get,
  568.  &random128_get_double};
  569.  
  570. static const gsl_rng_type random256_libc5_type =
  571. {"random256-libc5",            /* name */
  572.  0x7fffffffUL,            /* RAND_MAX */
  573.  0,                /* RAND_MIN */
  574.  sizeof (random256_state_t),
  575.  &random256_libc5_set,
  576.  &random256_get,
  577.  &random256_get_double};
  578.  
  579. static const gsl_rng_type random_bsd_type =
  580. {"random-bsd",            /* name */
  581.  0x7fffffffUL,            /* RAND_MAX */
  582.  0,                /* RAND_MIN */
  583.  sizeof (random128_state_t),
  584.  &random128_bsd_set,
  585.  &random128_get,
  586.  &random128_get_double};
  587.  
  588. static const gsl_rng_type random8_bsd_type =
  589. {"random8-bsd",            /* name */
  590.  0x7fffffffUL,            /* RAND_MAX */
  591.  0,                /* RAND_MIN */
  592.  sizeof (random8_state_t),
  593.  &random8_bsd_set,
  594.  &random8_get,
  595.  &random8_get_double};
  596.  
  597. static const gsl_rng_type random32_bsd_type =
  598. {"random32-bsd",            /* name */
  599.  0x7fffffffUL,            /* RAND_MAX */
  600.  0,                /* RAND_MIN */
  601.  sizeof (random32_state_t),
  602.  &random32_bsd_set,
  603.  &random32_get,
  604.  &random32_get_double};
  605.  
  606. static const gsl_rng_type random64_bsd_type =
  607. {"random64-bsd",            /* name */
  608.  0x7fffffffUL,            /* RAND_MAX */
  609.  0,                /* RAND_MIN */
  610.  sizeof (random64_state_t),
  611.  &random64_bsd_set,
  612.  &random64_get,
  613.  &random64_get_double};
  614.  
  615. static const gsl_rng_type random128_bsd_type =
  616. {"random128-bsd",        /* name */
  617.  0x7fffffffUL,            /* RAND_MAX */
  618.  0,                /* RAND_MIN */
  619.  sizeof (random128_state_t),
  620.  &random128_bsd_set,
  621.  &random128_get,
  622.  &random128_get_double};
  623.  
  624. static const gsl_rng_type random256_bsd_type =
  625. {"random256-bsd",        /* name */
  626.  0x7fffffffUL,            /* RAND_MAX */
  627.  0,                /* RAND_MIN */
  628.  sizeof (random256_state_t),
  629.  &random256_bsd_set,
  630.  &random256_get,
  631.  &random256_get_double};
  632.  
  633. const gsl_rng_type *gsl_rng_random_libc5    = &random_libc5_type;
  634. const gsl_rng_type *gsl_rng_random8_libc5   = &random8_libc5_type;
  635. const gsl_rng_type *gsl_rng_random32_libc5  = &random32_libc5_type;
  636. const gsl_rng_type *gsl_rng_random64_libc5  = &random64_libc5_type;
  637. const gsl_rng_type *gsl_rng_random128_libc5 = &random128_libc5_type;
  638. const gsl_rng_type *gsl_rng_random256_libc5 = &random256_libc5_type;
  639.  
  640. const gsl_rng_type *gsl_rng_random_glibc2    = &random_glibc2_type;
  641. const gsl_rng_type *gsl_rng_random8_glibc2   = &random8_glibc2_type;
  642. const gsl_rng_type *gsl_rng_random32_glibc2  = &random32_glibc2_type;
  643. const gsl_rng_type *gsl_rng_random64_glibc2  = &random64_glibc2_type;
  644. const gsl_rng_type *gsl_rng_random128_glibc2 = &random128_glibc2_type;
  645. const gsl_rng_type *gsl_rng_random256_glibc2 = &random256_glibc2_type;
  646.  
  647. const gsl_rng_type *gsl_rng_random_bsd    = &random_bsd_type;
  648. const gsl_rng_type *gsl_rng_random8_bsd   = &random8_bsd_type;
  649. const gsl_rng_type *gsl_rng_random32_bsd  = &random32_bsd_type;
  650. const gsl_rng_type *gsl_rng_random64_bsd  = &random64_bsd_type;
  651. const gsl_rng_type *gsl_rng_random128_bsd = &random128_bsd_type;
  652. const gsl_rng_type *gsl_rng_random256_bsd = &random256_bsd_type;
  653.  
  654.  
  655.  
  656.  
  657.